1 Relations across an oblique shock

Reference: P620-624 in John Anderson, “Fundamentals of Aerodynamics”,McGraw-Hill Education , 2017.

Control equations across the oblique shock are as follows,

\[ \begin{aligned} \rho_{1} u_{1} &= \rho_{2} u_{2} \\ w_{1} &= w_{2} \\ p_{1} + \rho_{1} u^2_{1} &= p_{2} + \rho_{2} u^2_{2} \\ h_{1} + \frac{u^2_{1}}{2} &= h_{2} + \frac{u^2_{2}}{2} \\ \end{aligned} \]

Geometry relations are as follows,

\[ \begin{aligned} M_{n,1} &= M_{1} \sin{\beta} \\ M_{n, 2}^{2} &= \frac{1+[(\gamma-1) / 2] M_{n, 1}^{2}}{\gamma M_{n, 1}^{2}-(\gamma-1) / 2} \\ \frac{\rho_{2}}{\rho_{1}} &=\frac{(\gamma+1) M_{n, 1}^{2}}{2+(\gamma-1) M_{n, 1}^{2}} \\ \frac{p_{2}}{p_{1}} &= 1+\frac{2 \gamma}{\gamma+1}\left(M_{n, 1}^{2}-1\right) \\ \frac{T_{2}}{T_{1}} &= \frac{p_{2}}{p_{1}} \frac{\rho_{1}}{\rho_{2}} \\ \tan{\theta} &= 2 \cot{(\beta)} \frac{M_{1}^{2} \sin^{2}{(\beta)}-1}{M_{1}^{2}(\gamma+\cos{(2 \beta)})+2} \\ M_{2} &= \frac{M_{n, 2}}{\sin{(\beta-\theta)}} \\ \end{aligned} \]

Python script of \(\theta-\beta-M\) relation is as follows,

def theta_beta_M(beta,theta,M1):
    from numpy import sin,cos,tan
    return (M1**2 * sin(beta)**2 - 1) / (M1**2*(GAMMA+cos(2*beta))+2) * 2/tan(beta) - tan(theta)

Use root-finding to find \(\beta\) if \(\theta\) is given,

from scipy.optimize import root
res = root(theta_beta_M,ramp_angle,args=(ramp_angle,Ma_infty))
beta = res.x
print("Beta is %.14f"%(beta/np.pi*180.0))

A python script to handle both the ramp case and the incident shock case is as follows,

import numpy as np

GAMMA=1.4
R=287.058

def theta_beta_M(beta,theta,M1):
    GAMMA=1.4
    from numpy import sin,cos,tan
    return (M1**2 * sin(beta)**2 - 1) / (M1**2*(GAMMA+cos(2*beta))+2) * 2/tan(beta) - tan(theta)

class ShockRelation:

    def __init__(self,M1,rho1,p1):
        from numpy import sqrt
        self.M1 = M1
        self.rho1 = rho1
        self.p1 = p1
        self.T1 = self.p1/(self.rho1*R)
        self.c1 = sqrt(GAMMA*R*self.T1)
        self.u1 = self.M1 * self.c1
        self.v1 = 0.0
        print("Velocities before the shock are [%23.15e, %23.15e]"%(self.u1,self.v1))
        self.beta = 0.0
        self.theta = 0.0
        self.Mn1 = 0.0
        self.p_ratio = 0.0
        self.rho_ratio = 0.0
        self.T_ratio = 0.0
        self.p2 = 0.0
        self.rho2 = 0.0
        self.T2 = 0.0
        self.c2 = 0.0
        self.Mn2 = 0.0
        self.M2 = 0.0
        self.u2 = 0.0
        self.v2 = 0.0

    def setBeta(self,beta):
        self.beta = beta

    def setTheta(self,theta):
        self.theta = theta

    def setBetaByTheta(self):
        if (self.theta==0.0):
            exit("self.theta(==0.0) is not set!")
        from scipy.optimize import root
        res = root(theta_beta_M,self.theta,args=(self.theta,self.M1))
        self.beta = res.x
        print("Beta is %23.15e"%(self.beta/np.pi*180.0))

    def setThetaByBeta(self):
        if (self.beta==0.0):
            exit("self.beta(==0.0) is not set!")
        from numpy import arctan
        self.theta=arctan(theta_beta_M(self.beta,0.0,self.M1))
        print("Theta is %23.15e"%(self.theta/np.pi*180.0))

    def setShockPostState(self):
        if (self.beta==0.0):
            exit("self.beta(==0.0) is not set!")
        from numpy import sin,cos,sqrt
        self.Mn1 = self.M1 * sin(self.beta)
        self.p_ratio = 1.0 + 2.0*GAMMA/(GAMMA+1.0) * (self.Mn1**2-1.0)
        self.p2 = self.p_ratio*self.p1
        print("P2/P1 = %23.15e, P2 = %23.15e"%(self.p_ratio,self.p2))
        self.rho_ratio = ( (GAMMA+1.0) * self.Mn1**2 ) / ( 2.0 + (GAMMA-1.0) * self.Mn1**2 )
        self.rho2 = self.rho1 * self.rho_ratio
        print("Rho2/Rho1 = %23.15e, Rho2 = %23.15e"%(self.rho_ratio,self.rho2))
        self.T_ratio = self.p_ratio / self.rho_ratio
        self.T2 = self.T1 * self.T_ratio
        print("T2/T1= %23.15e, T2 = %23.15e"%(self.T_ratio,self.T2))
        self.Mn2 = ( 1.0 + ((GAMMA-1.0)/2.0)*self.Mn1**2 ) / ( GAMMA*self.Mn1**2 - (GAMMA-1.0)/2.0 )
        self.M2 = self.Mn2 / sin(self.beta-self.theta)
        print("Ma behind the shock is %23.15e"%self.M2)
        print("Normal Ma across the shock are [%.3f,%.3f]"%(self.Mn1,self.Mn2))
        self.c2 = sqrt(GAMMA*R*self.T2)
        vel2 = self.M2 * self.c2
        self.u2 = vel2 * cos(self.theta)
        self.v2 = vel2 * sin(self.theta)
        print("Velocities behind the shock are [%23.15e, %23.15e]"%(self.u2,self.v2))

if __name__ == '__main__':
    # Examples
    ## Ramp case
    Ma_infty = 2.0
    rho_infty = 1.17683
    p_infty = 101325.0
    shock_relation = ShockRelation(M1=Ma_infty,rho1=rho_infty,p1=p_infty)
    ramp_angle = 15.0/180.0*np.pi
    shock_relation.setTheta(ramp_angle)
    shock_relation.setBetaByTheta()
    shock_relation.setShockPostState()

    ## Incident shock case
    Ma_infty = 2.15
    rho_infty = 2.445716025859303e-03
    T_infty = 288.15
    p_infty = rho_infty * R * T_infty
    shock_relation = ShockRelation(M1=Ma_infty,rho1=rho_infty,p1=p_infty)
    beta = 30.8/180.0*np.pi
    shock_relation.setBeta(beta)
    shock_relation.setThetaByBeta()
    shock_relation.setShockPostState()

2 Rankine-Hugoniot relation

Use RH relation to find the shock speed. Rankine-Hugoniot condition is based on the same control equations. Check Rankine–Hugoniot conditions: Shock_condition in Wikipedia.

\[ u_{s}=u_{1}+c_{1} \sqrt{1+\frac{\gamma+1}{2 \gamma}\left(\frac{p_{2}}{p_{1}}-1\right)} \]

where \(c_{1}=\sqrt{\gamma p_{1}/\rho_{1}}\) is the speed of sound in the fluid at upstream conditions

Shock relations

2020-09-24

1 Relations across an oblique shock

2 Rankine-Hugoniot relation